Scale-Free topologies and Activatory-Inhibitory interactions 
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A simple model of activatory-inhibitory interactions controlling the activity of agents (substrates) 
through a "saturated response" dynamical rule in a scale-free network is thoroughly studied. Af- 
ter discussing the most remarkable dynamical features of the model, namely fragmentation and 
multistability, we present a characterization of the temporal (periodic and chaotic) fluctuations of 
the quasi-stasis asymptotic states of network activity. The double (both structural and dynamical) 
source of entangled complexity of the system temporal fluctuations, as an important partial aspect 
of the Correlation Structure- Function problem, is further discussed to the light of the numerical 
results, with a view on potential applications of these general results. 

PACS numbers: 89.75.Fb, 05.45.-a 



o : 

x> : 
i 

Si 

> ■ 
on : 

(N . 
O ■ 
■ 

o : 
o ■ 



i , 



X 



Many real networks are complex and heteroge- 
neous. In this paper, we study a dynamics that 
generically describes biological processes that 
take place on complex architectures as metabolic 
reactions and gene expression. We capitalize on 
the theory of nonlinear dynamical systems to un- 
cover the topological and dynamical patterns of a 
Michaelis-Menten like dynamics coupled to a net- 
work that is complex, directed and highly skewed. 
The results indicate that such patterns can exist 
in the form of periodic and chaotic orbits reveal- 
ing interesting properties at both local and global 
levels. Moreover, the dynamics on top of the 
substrate networks yields topologically complex 
substructures (islands or clusters) whose struc- 
tural charateristics are analysed. This analysis of- 
fers interesting results on the interplay Structure- 
Function. We round off our study by discussing 
possible implications of heterogeneous topologies 
on biological processes at the cellular level. 



I. INTRODUCTION 

Nonlinear lattices, i.e. spatially discrete many-body 
systems with nonlinear interactions, are currently the 
subject of a considerable multidisciplinary interest, not 
only among a wide variety of Physics subdisciplines and 
technologies Q, but also in Biomolecular Chemistry, 
Cell physiology, Theoretical Biology, Social Sciences and 
other fields Q|. There is a common basic interest in the 
understanding of the many aspects of the correlation be- 
tween Structure and Function in systems made up of dis- 
cretely many nonlinearly interacting components. 

While in Physics applications the interactions among 
constituents (atoms, magnetic moments, . . . ) usually de- 
pends on the geometrical distance between their posi- 



tions in real space, in applications outside "fundamental 
physics" the space of interactions is abstracted, so that 
proximity between components (say i and j) is measured 
as length of the path of interactions given by a connec- 
tivity matrix dj. In other words, the graph (or network) 
of interactions, and not anymore the real space, fixes the 
relevant geometry related to the function (dynamics) of 
the system. This does not preclude the applicability of 
Statistical and Field theory methods to the study of non- 
linear lattices outside the traditional physics subdisci- 
plines; on the contrary, a proper use of these approaches 
often throws considerable light on some important issues 
currently addressed Q. 

Though lattice disorder effects on nonlinear dynam- 
ics of macroscopic systems have their own tradition, the 
most usual case in physics is that of homogeneous (ei- 
ther pure random or regular) networks. However, re- 
cent confluent studies on the Structure of interactions in 
a large variety of technological (communication, power 
grid) as well as biomolecular (protein-protein interaction, 
gene regulation, cell metabolism), ecological (trophic 
networks, mutualism) and socio-economic systems have 
shown the overabundance of highly inhomogeneous struc- 
tures among "real world" interacting systems. 

Homogeneity of the interactions structure means that 
almost all nodes are topologically equivalent, like in a 
regular lattice or in a random Erdos-Renyi graph, thus 
showing a density distribution function of the degree of 
connectivity localized around a mean value with a well- 
defined average of quadratic fluctuations. If k denotes 
the degree (number of interactions of a given node) and 
P(fc) denotes its density distribution function, an inho- 
mogeneous network shows for P(k) a power law (often 
truncated). The absence of a characteristic scale in the 
connectivity patterns (scale-free networks) manifests it- 
self in the presence of a few number of nodes (named 
hubs) connected to very many nodes, and a larger num- 
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ber of poorly connected nodes. The complex character of 
the structure of the interactions couples to the dynamical 
complexity which emerges from the nonlinear character 
of the interactions, so that generally one may say that the 
Structure- Function correlation problem in real networks 
has at least two sources of entangled complexity. 

The model that we analyze is introduced in section 
ITT1 It tries to capture general ingredients of this entan- 
gled complexity in a relevant kind of nonlinear dynamics: 
Activation/Inhibition (AI) competing interactions with 
a "saturated response" rule for the rate of activation 
(see Figure 1). This kind of dynamics is often called 
Michaelis-Menten Q, Holling 0|, or Langmuir Q rule. 
The interacting units sit on a lattice which is both small- 
world (i.e. short mean path length) and scale-free. For 
this we use the Barabasi- Albert |S( network. Afterwords, 
some basic general features of the model are discussed, 
namely the network fragmentation in subclusters (or is- 
lands) of collective dynamics (jll Afl . and the generic types 
of asymptotic behaviours coexisting in the phase space of 
collective dynamics as well as the observed bifurcations in 
phase portrait upon parameter variations (|IIB|I . These 
basic consequences of the AI competition on the complex 
network are prevalent for values of the ratio AI ranging 
from 1 to 6. Finally, in section lll Cl the bifurcations found 
are explained in terms of the Floquet analysis of the so- 
lutions. 

Section IlIII mainly reports on the statistical character- 
ization of both the dynamical behaviours observed (sec- 
tion llll A"|) and the structural characterization of the dy- 

We perform an extensive 



namical islands fsection flll Bll 
exploration of the parameter space, employing different 
initial conditions and substrate network realizations, in 
order to find the conditions for the existence of chaotic 
and periodic behavior as well as to fully characterize the 
main topological characteristics of the dynamical islands. 
Finally, in section lTlI Cl we identify those substructures of 
the dynamical islands that are relevant for the dynamical 
evolution of the system. 

The concluding section llVl summarizes the main con- 
clusions of our work, along with some prospective re- 
marks on likely applications of the model, and the poten- 
tial use of these techniques in the study of particularly 
interesting real-world biological networks. 



II. THE MODEL. BASIC DYNAMICAL 
FEATURES. 

As stated in the introductory section ^ we introduce 
here a model of Activatory/Inhibitory interactions regu- 
lating the activity gi{t) (i = 1, .., AT), of N constituents 
(e.g. agents, substrates), with N generally being a large 
number. The real functions of time gi(t) are each one 
attached to a node of a graph with adjacency matrix CV,- 
[N x N). The matrix element is non-zero, Cy ^ 0, only 
if the rate of variation of the i-th node activity, gi(t), 
depends on the activity gj of the j-th node (interaction 



i *— j). Different realizations of the Cy matrix are con- 
structed using the method of Barabasi and Albert 0, in 
order to ensure two seemingly universal characteristics of 
many recently studied networks in Biological and Social 
sciences and other fields 



(a) Small-world, meaning that the mean distance (min- 
imal length of the interaction path), (%), between 
pairs of nodes goes at most as log N, for large val- 
ues of N. 

(b) Scale-free, meaning that the density distribution 
function P(k) of the degree (connectivity) of nodes 
scales as P(k) ~ A: -7 , with 7 = 3. Other values of 
7 (2 < 7 < 3) were also analyzed by using suitably 
tested modifications of the Barabasi- Albert prefer- 
ential attachment rule [Io|. 

The interaction (i <— j) can be either activatory (exci- 
tatory) or inhibitory; correspondingly we define the inter- 
action matrix element Wij to be +1 or —1, respectively 
(and Wij — whenever Cy = 0), and call p the fraction, 
among non-zero elements, of negative signs (note that 
while C'ij is a symmetric matrix, Wij is not in general). 
Moreover, the sign distribution of elements is taken uni- 
form in the set of (approx. (k)N/2) links of the network 
realization. 
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FIG. 1: Saturated functional response (h — 4). 

The dynamics of the nodes activity vector G(t) — 
{gi(t)} (with i — 1, . . ., N) that we consider is such that 
only in the pressence of excitatory neighbours activity 
could gi possibly be non null, otherwise gi decays to zero 
with an exponential rate: 



dGjt) 
dt 



-G(t) + aF(WG(t)), 



(1) 



where F is a nonlinear vector function whose argument 
is the product of the interaction matrix W and the ac- 
tivity vector G, and a (> 0) accounts for the strength of 
the interaction. The function F is a saturated response 
function (see figure defined as: 



FIG. 2: (Color online) Two examples of network fragmentation. The nodes of the networks are clasified in: (i) dynamical nodes 
(red), (it) stationary nodes with nonzero activity (blue), (in) stationary nodes with zero activity beloging to T>q ©X>i (yellow) 
and (iv) remaining nodes with zero activity (white). Note that the white central nodes in (b) act as the frontier between the 
dynamical island and the steady nonzero activity one. 



where $(x) is a function denned as $(x) = x if x > 
and zero otherwise. In our numerical studies of the 
model we have fixed the value of the parameter a = 3, 
and varied the parameters < p < 1 and < h < 10. 
One can easily realize that the solutions for non-negative 
initial conditions remain bounded for all t > 0. Indeed, 
as the nonlinear term in is bounded above by a, one 
obtains that c/i < whenever gi > a. Also, if gi = 
then Fj(WG) > 0, so that the activities cannot become 
negative. 

The above dynamics can be regarded, e. g., as a gen- 
eralization of some simplified and coarse-grained genetic 
models, referred to as Random Boolean Networks [TT| . 
where Boolean rules are implemented. These have been 
extensively used to study networks at various levels of bi- 
ological organization, and have provided useful insights 
later supported by experiments [ll], • Equation in- 
corporates the experimental observation of a continuous 
range of activity levels While linear models have 

been successful for the reconstruction of the interaction 
networks from experimental data |l4j . nonlinear models 
like Eqs. are expected to be more appropriate for a 
quantitative description of the dynamics. 

The dynamics JIJ of a two-agent (dimer) model has 
been considered in reference [l3| , in the context of virus- 
cell interactions in bacteria and general gene regulatory 
activity models, where a rich repertoire of behaviours, 
like multi-stability (multiple attractors in phase space) 
was reported. A preliminary study of the behaviour of 
(ffl) on small-world scale-free networks can be found in 
fl5| , where the interested reader can find a more detailed 
account of the numerical techniques used in the char- 
acterization of the different dynamical regimes. In this 
section, we review some remarkable general features of 



the network dynamics. 



A. Activation and Inhibition interplay. 
Fragmentation. 

A brief look at equation easily reveals that for any 
value of the parameters p and h the state of inactivity, 
G = 0, is always a solution. As a matter of fact, for 
h = 0, or h 7^ but p = 1, this is the unique asymptotic 
solution (global attractor in the phase space) for all pos- 
sible non-negative initial conditions. However, for h =/= 
andp 7^ 1 other asymptotic solutions, with islands of pos- 
itive activity, generically coexist with the rest state. The 
term islands denotes here subnetworks that are intercon- 
nected through nodes which have evolved to null activity, 
so that their dynamics are effectively disconnected. 

The fragmentation of the network dynamics into dis- 
connected islands is a generic feature of AI interactions, 
as the following considerations suggest. Let us call V the 
set of nodes whose activities, for a given initial condition 
G(t = 0), asymptotically vanish. It is easy to see that, 
irrespective of the initial condition, this set is generically 
non-empty. 

Indeed, if a node i is such that Wij = — 1 or for all j, 
then its activity gi(t) will tend to zero. Let us call £>o the 
set of these nodes, and note that its measure Q2k P(k)p k ) 
is a non-zero increasing function of p. Now, call 2?i the 
set of nodes I such that their positive Wij occur for j's 
in T>q, and so on ... Due to the small- world property, 
there are in fact very few relevant V n (n = 0, 1, . . .) sets. 
Its union T>* = (J T> n is easily seen to have a non-zero 
measure which increases with p. 

The nodes of T>* are structurally {i.e. irrespective of 
initial conditions) inactive. Depending on the initial con- 
dition, the set T> may include other nodes not contained 
in T>* , namely those nodes that evolve to inactivity due to 
the initial condition (dynamically, instead of structurally, 
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inactive): See e.g. the white nodes in figure [2 where we 
show two small networks of N = 50 nodes to allow a sim- 
ple visualization of the sets T>* and T>. In other words, 
the measure of T> may in general be (much) larger than 
the measure of the "structurally dead" nodes V* . 

From the previous considerations, wether or not the set 
T> percolates the network realization, leaving out islands 
of disconnected activity, is an event that clearly depends 
on both the parameter p and the initial conditions. But 
also the discussion correctly suggests that fragmentation 
of the network into subclusters with independent tempo- 
ral evolution is a generic (non-zero measure) feature. Our 
numerics, which are extensive in the sense of (both, net- 
work realizations and initial conditions) large sampling, 
convincingly corroborate this assertion. Figures and 01 
show two islands of periodic and chaotic activity, respec- 
tively, as well as the temporal evolution of g%(t) for some 
of their constituent's nodes (see the next section for a 
more detailed discussion of the figures). 



B. Temporal fluctuations of asymptotic solutions. 

The asymptotic dynamics of equation (JJJ was studied 
in |15|. We summarize here the most salient features of 
the phase portrait of the collective dynamics. 

The presence of inhibitory interactions makes possible 
the existence of instabilities in the fixed point solutions 
{i.e. states of constant activities, gi(t) — g*, let us say 
chemostasis regime) of evolution equations Q. Using 
linear stability analysis techniques, these "typical" insta- 
bilities are characterized as Hopf bifurcations (either di- 
rect or often inverse), where attractors of exactly periodic 
collective activities, <ft(i) = gt(t + T), are born out from 
the unstable fixed points. The inverse period (frequency) 
u> = 1 /T of a periodic attractor changes with parameter 
and is naturally dependent on each specific island real- 
ization. A sampling over different initial conditions and 
network realizations shows that lo is smoothly peak- dis- 
tributed around a value of order unity (the characteristic 
time scale of activity decay in the absence of interactions) 
with a decaying queue slightly biased to higher frequency 
values [l5| . 

One easily observes that these periodic attractors, in 
turn, typically experience also period doubling instabil- 
ities, and through the well-known universal scenario of 
(successive) period doubling bifurcation cascade, the on- 
set of chaotic attractors takes place in the phase portrait 
of the network dynamics. To help visualization of the 
generic types of asymptotic network dynamics, we rep- 
resent in figure [5] the bifurcation diagram for a typical 
attractor. At different values of the (Michaelis-Mcnten) 
parameter h, and constant values of a = 3,p — 0.7, we 
plot the activity of an individual node at the instant when 
its time derivative vanishes. Thus, a single branch in the 
figure indicates stationary activity, two branches indicate 
a periodic attractor, etc. We also plot in figure the 
largest Lyapunov exponent A on the attractor, so to al- 



low discerning between chaotic (positive A) and eventual 
regular quasiperiodic evolutions (A = 0). 

A similar bifurcation diagram for a different network 
realization is shown in figure where one can appreci- 
ate (see inset) a commonly found bifurcation (though it 
appears much less often than period doubling), namely 
period tripling bifurcation. Its characterization will be 
made below in the next subsection where the Floquct 
analysis of periodic attractors is presented. 

To illustrate the aspect of typical periodic fluctuations 
we show in figure |31 some examples of the temporal ac- 
tivity gi(t) of different nodes inside an island of syncro- 
nised activity from a representative system. Note that 
the abundance of out of phase oscillations of neighbours 
activity is a natural consequence of inhibitory interac- 
tions. Horizontal lines in insets indicate the average level 
cji of node activity. We see that in some of the island 
nodes the amplitude of the oscillation is small compared 
to <ji (see e.g. top rightmost and bottom leftmost in- 
sets); while in others they are of comparable size, even 
to the point that lowest levels of activity can reach a null 
value, before activity is triggered again after inhibiting 
neighbors activity decreases enough. 

An analogous visualization of chaotic temporal fluctu- 
ations of the activities in a cluster is shown in figure 01 
Here again we see nodes {e.g. top left inset) where the 
size of activity fluctuations is less than 1 per cent of the 
average level g~i. Most remarkable, there are nodes (like 
the one in bottom left inset) which remain inactive most 
of the time intermittently experiencing spikes of short 
duration activity. This amazing variability of individ- 
ual node temporal activity on the chaotic attractors is a 
generic feature of the network dynamics. The existence of 
spike behaviour of individual nodes activity suggests cor- 
rectly that eventual variations of parameters like h may 
lead to permanent inactivity of some particular nodes, 
so providing a straightaway decreasing of the dynamical 
cluster size or, the other way around, the activation of 
inactive nodes in the frontier. 

It is important to note that, for a fixed set of parameter 
values and a given network realization, there are gener- 
ally several different attractors coexisting in the phase 
space portrait of the network dynamics, each one having 
its own basin (of attraction) of initial conditions. Multi- 
stability appears as a generic consequence of the exci- 
tatory/inhibitory interplay. Importantly also, there can 
be very many unstable periodic trajectories (often entan- 
gled) flowing in between basins of attractions. The exci- 
tatory/inhibitory competition is also responsible for the 
appearance of temporally complex (positive Lyapunov 
exponent) aperiodic evolutions, associated to the bifurca- 
tion cascade scenario. As we will show in section llH Al thc 
manifestation of fluctuating (either periodic or chaotic) 
temporal behaviours takes importance when inhibitory 
links predominate, though not too much, over excitatory 
ones. 
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C. Floquet analysis of the periodic attractors. 

As shown in the bifurcation diagrams of figures [S] and 
periodic solutions of the network dynamics often be- 
come unstable under variations of the model parameters. 
In order to characterize these instabilities in a precise 
manner, one may perform the linear stability analysis of 
the periodic orbits (see, e.g. 01) near the bifurcation 
points. 

For this we consider small perturbations of the dy- 
namical variables, 5g(to) — {Sgi(to)}, and compute their 
evolution over the period T of the periodic orbit. The 
evolution of these small perturbations (vectors in tan- 
gent space) follows the (linear) dynamics obtained by lin- 
earizing equation JQ) around the periodic orbit {(ji{t)} = 
{gi(t + T)},i.e., 



dSgjt) 
dt 



a ■ ASg 



where the matrix A is obtained as 



■ w, 



h3 



(3) 



(4) 



a+h-^[j: k w hk g k ]r 

and 0[x] denotes the (Heaviside) step function. Note 



that the above equation is only valid when the sum of 
the inputs (activatories and inhibitorics) which receives a 
node from its neighbours is nonzero. Hence, the Floquet 
analysis is performed for each dynamical cluster found 
and not for the whole network. 

The so-called Floquet (or monodromy) matrix T of the 
periodic solution {gt (t) } is defined as the linear operator 
in tangent space that maps the initial perturbation at to, 
5g(to), onto the perturbation at to + T 



5g{t + T) = T5g{to) 



(5) 



The Floquet matrix J- is obtained by numerical inte- 
gration of the linearized equation © over a period T for a 
basis of initial conditions in the tangent space. The spec- 
trum of eigenvalues of this matrix provides the informa- 
tion on the linear stability of the periodic solution. Note 
that because J- is & real matrix, if a Floquet eigenvalue fi 
is a complex number, then its complex conjugate fl also 
belongs to the Floquet spectrum. Also, because solutions 
of autonomous differential equations can be shifted in the 
time t direction, their Floquet matrix always has unity 
as an eigenvalue, say fi — 1, with associated eigenvector 
{p,-(to)}. The solution is linearly stable if all the other 
eigenvalues (i 1 = \fi 1 \exp(i9 J ) are in the interior of the 
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FIG. 5: Example of bifurcation diagram (island size: 14; N = 
60; p = 0.8). One can appreciate an inverse Hopf bifurcation 
and several (direct and inverse) period doubling bifurcation 
cascades. The maximum Lyapunov exponent A is plotted in 
the lower part. 




FIG. 6: Example of bifurcation diagram (island size: 12; N = 
60; p = 0.8) showing (see inset) a period tripling bifurcation. 
The maximum Lyapunov exponent A is plotted in the lower 
part. 
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FIG. 7: Floquet spectra, (a) Period doubling bifurcation 
in an island of size 14 at h = 1.57. (b) Naimark-S acker 
bifurcation (at rational Floquet angle 6 = ±27r/3) in an island 
of size 12 at h=2.44 (the same as used in the diagram of figure 
0. For both N = 60 and p = 0.8. 

unit circle of the complex plane, i.e. \(jp\ < 1 for j ^ 1. A 
periodic solution becomes unstable when a Floquet eigen- 
value (or a pair of complex conjugate eigenvalues) crosses 
the unit circle. The associated Floquet eigenvector indi- 
cates the direction in tangent space where perturbations 
will grow exponentially away from the solution. 

In figure Efa) we plot the Floquet spectrum of a peri- 
odic attractor at a period doubling bifurcation. As seen 
in the figure, a Floquet eigenvalue crosses the unit circle 
at — 1. In figureCJb) we plot the Floquet spectrum of the 
periodic attractor of figureHJlat h — 2.44, where the inset 
suggested that a period tripling bifurcation may occur. 
We see a complex conjugate pair of Floquet eigenvalues 
exiting the unit circle at angles 9 — ±2tt/3. In general, 
for generic irrational values of 8/n this type of bifurcation 
(called Naimark-Sacker or generalized Hopf bifurcation) 
gives rise to a quasiperiodic attractor whose trajectories 
fill densely a two-frequency torus. However, as a generic 
feature of our model, the two frequencies of the new at- 



tractor are in a commensurate ratio (2 : 3), so that the 
new stable trajectory has a period of 3T. 

In terms of how often different types of bifurcation oc- 
cur in the network dynamics, as inferred from our (non- 
exhaustive, but significant at the scales considered) sam- 
pling of initial conditions and network realizations, one 
may say that period doubling cascades and, less often, 
commensurate Naimark-Sacker bifurcations have been 
generically found by varying the Michaelis-Menten pa- 
rameter h. But, besides the formal characterization of 
the dynamical instabilities observed, the Floquet analy- 
sis allows also to give an answer on a more general ques- 
tion, namely how temporal instabilities correlate with 
networking connectivity characteristics. Are there char- 
acteristic features discernible in the structure of instabil- 
ities? This point will be discussed further below in the 
next subsection. 



III. STATISTICAL CHARACTERIZATION OF 
DYNAMICAL REGIMES AND ISLANDS 
STRUCTURE. 

As noted before, the dynamics of the system is deter- 
mined by only two parameters, h and p. The behaviour 
of the system described by equation Q on the underlying 
network is very rich and one can have steady, periodic or 
chaotic states as well as fragmentation. In this section, 
we analyze in more details the system's phase diagram 
as well as how the dynamical regimes couple to the lo- 
cal structural properties of the underlying network and 
dynamical islands. 



A. Density distribution functions of dynamical 
regimes. 

In figure|S| we have represented the probability, P c haos, 
of ending up in a chaotic regime as a function of p for a 
network of N = 300 nodes and h = 4. This probability is 
given by the fraction of the total number of realizations 
(typically 10 3 different initial conditions over different 
network realizations for each value of p and h were used) 
in which at least one chaotic dynamics is observed. The 
figure also shows the corresponding probability, P per , for 
periodic orbits. As figure Ufa) clearly shows, there is a 
threshold value p = pth beyond which the network dy- 
namics is not robust under variations of the initial values 
of the ijj's. For values of p above pth ~ 0.25(5), two 
randomly chosen initial conditions can lead the system 
to disparate asymptotic regimes. Besides, the size of the 
system affects the value of P c haos, but the onset —and 
the end— of the chaotic phase seems to be N indepen- 
dent [13. 

Moreover, figure Ufa) constitutes a quantitative illus- 
tration of how the prevalence of fluctuating asymptotic 
regimes over chemo-stasis ones depends on the model 
parameter p. The sum of both functions, P per (p) + 
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FIG. 8: (Color online) (a) Probability, P c haos (Pper), that the system evolves to a chaotic (periodic) regime as a function of 
the probability of inhibitory interactions, p, for h — 4 and N = 300. (b) Phase diagram in the (p, ft)-parameter space of the 
chaotic dynamics of the system. Color code indicate the values of P c ha OS (N = 300) . 
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Pchaos{p)i gives the probability that the asymptotic state 
shows temporal variations of the activity vector (either 
regular or chaotic) as a function of p. These results give 
that in the range of values 0.5 < p < 0.8 regimes of tem- 
poral fluctuations occur more often than constant activ- 
ity regimes. This measure is maximized by values around 
p ~ 2/3 and, quite naturally, it increases with the value of 
the Michaelis-Menten parameter h, i.e. the slope at the 
origin of the saturated response function (see figure [T]l. 
Note that even larger values of p means overabundance 
of inhibitory interactions, which leads to the predomi- 
nance of the asymptotic rest state, while smaller values 
of p favour chemostatic equilibria. 
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FIG. 9: (a) Probability that a connected cluster of nodes 
displaying either chaotic or periodic behavior has a given size 
(in number of nodes forming the cluster), (b) Scaling of the 
mean cluster size with N. The parameters have been set to 
h = 4 and p = 0.7. 



The quantities P c haos and p t h depend on h. As we 
move to larger values of h, the strength of the inter- 
actions increases and hence it is expected that slight 
perturbations produce a behavior in which the fraction 
of nodes whose dynamical patterns are easily disturbed 
grows. This is indeed the illustrated in figure 

|H(b). The color-coded figure shows that as h is increased, 
the probability of having a chaotic phase grows, and that 
the onset of such chaotic patterns shifts to smaller values 
of p. This drift of pth is however bounded. For small 
enough values of p (even for very large h), most of the 
elements activate each other (Wij = 1 for a large fraction 
of pairs ij and ji) and hence the resulting dynamics is 
steady. In other words, the onset of chaotic regimes is al- 
ways located at a nonzero value of pth (the same applies 
to the right (decaying) part of P c h{p)> but in this case 
the activity falls down to zero). 
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FIG. 10: Probability that a node belonging to a dynamical 
island interacts with k other nodes of the island. Parameters 
were set to h — 4 and p — 0.7. 
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FIG. 11: Average clustering coefficient (C) as a function of the 
network size for the BA original network and the dynamical 
cluster. Note that while (C) in the BA network continuously 
decreases, for the dynamical island it saturates. See more 
details in the text. The results have been obtained using 
h = 4 and p = 0.7. 



dynamics appears as one of the most characteristic fea- 
tures of the model. 

Another interesting aspect is the elucidation of how 
the topological properties of the islands correlate with 
those of the underlying (original) network. To this end, 
we have further scrutinized the structure of the clusters 
and measured two topological quantities of interest. Fig- 
ure EI| shows the degree distribution of nodes belonging 
to dynamical islands for several system sizes. This prop- 
erty can be regarded as a global one and indicates that 
within the islands, the probability that a node has k links 
also follows a power-law, though with a more pronounced 
cut-off and a (slightly) different value for the exponent 7. 
More striking is the result depicted in figure^J where the 
average clustering coefficient (C) of the substrate (orig- 
inal) network and of the islands is plotted as a function 
of N. While for the BA network the clustering is vanish- 
ing as the network size grows, it seems that for dynam- 
ical islands its value saturates. This is quite interesting 
because, on the other hand, the value of the clustering 
coefficient is very large and comparable to measures of 
real systems where the kind of dynamics explored here 
applies, for instance, biological networks |T(| . 

That is, the structure of dynamical islands correctly re- 
produces several of the most important topological fea- 
tures observed in biological networks and not captured 
by current network models. Namely, the heterogeneous 
distribution of connections, a high average clustering and 
the independence of (C) with respect to the system size. 
This result points to the conjecture that several topolog- 
ical properties observed in systems driven by AI interac- 
tions where nodes are themselves (nonlinear) dynamical 
units may be biassed by their own dynamics. In other 
words, what we actually see is the result of the activity 
showed up by a smaller "dynamical" network whose local 
topological properties greatly differ from those of a larger 
substrate network that we don't "see" because many of 
its components are simply off. This, in fact, may be the 
case of biological systems where structure and dynamics 
are indissoluble linked 0, 0] . 



B. Dynamical island structure 



C. Structure inside Dynamical Islands 



We next focus on the topological characterization of is- 
lands of dynamical units. We first analyze how the clus- 
ter size distribution of islands of nodes displaying either 
periodic or chaotic activity scales with the system size. 
Figure Eta) represents the probability that an island has 
a given size for several networks made up of a number 
of nodes ranging from 50 to 800. Clearly, the size distri- 
bution shows an average value that changes as N grows. 
A closer look at the figure (see figure Efb)) reveals that 
the mean cluster size scales with N and that about 17% 
of the nodes, in average, exhibits nonzero activity. This 
confirms what we have discussed in section III AI about 
the measures of the sets V* and V, namely, that the 
fragmentation of the network into islands of independent 



The above findings on new (dynamically) emergent 
characteristics of the islands structure motivates the 
question of wheter these clusters have an internal orga- 
nization or hierarchy among its constituents. It is widely 
known that when one deals with problems where the net- 
work topology (scale-free) is the only degree of complex- 
ity of the problem the answer to this question is usu- 
ally based on the presence of highly conected nodes (the 
hubs). This is the case when linear evolution equations 
are studied on top of complex networks like epidemic or 
rumour spreading, traffic and communication problems 
[l8t fl9t l20j. However, our case is not so simple and the 
nonlinear excitatory/inhibitory dynamics between the el- 
ements of the network plays a crucial role in determining 
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FIG. 12: (Color online) (a) and (d) show the components of the vector (Sg ) (see text) for two dynamical islands at the 
critical point of a Naimark-S acker bifurcation at h = 2.44 (a) and a period doubling one at h = 2.629 (d). (b) and (e) show 
the distribution (green region) of the nodes with non null component of (Sg ) in (a) and (d) respectively. Finally, (c) and 
(f ) show the evolution of the Floquet eigenvalue fi* as a function of the forcing amplitude A applied to different nodes of the 
dynamical islands (b) and (e) respectively. For both islands the susbstrate network was of N = 60 nodes whith a fraction of 
80% of inhibitory interactions (p = 0.8). 



which nodes are governing the evolution of the system. 
Moreover, the high clustering found for the dynamical 
clusters points out that this leading role is not played 
by isolated nodes but by small substructures inside the 
dynamical islands. This concept is not new, the problem 
of finding small relevant substructures inside large net- 
works, usually called "motifs" j^, has been studied in 
different ways in the field of biological networks. 

It is indeed very revealing to pay attention to the net- 
worked structure of the unstable manifold, which is given 
in the linear regime of small perturbations by the Flo- 
quet unstable eigenvectors. For this purpose, we look at 
the behavior of the components of the dynamical islands 
when a bifurcation (either period doubling or Naimark- 
Sacker type) occurs. In these critical points, it is possible 
to get a deeper insight into what is going on in the dy- 
namical islands I by looking at the Floquet eigenvector 
responsible for the bifurcation, Sg (to) — {5g*(t )}, cor- 
responding to the Floquet eigenvalue which reaches the 
unit circle. In particular, integrating equation © for the 
initial condition 5g (to) we can compute the following 
vector 

(S?) = {{Sg*}} = ^ £ \Sgt\dt^ . (6) 

The components of this vector measure, for each node, 
the average (over a period T of the old solution) distance 
of the new solution after the bifurcation point from the 
old periodic solution. Note that a zero component of this 



vector at a node k, means orthogonality of the single-site 
perturbations at that node with respect to the unsta- 
ble direction in tangent space. In other words, by look- 
ing at the components of the vector © we can identify 
those nodes that are more affected by the perturbation 
that leads the system to instability. In figures IT^T a) and 
I12f d) we show this quantity for the same two islands 
(relatively small, but still representative) whose Floquet 
spectra are given in Fig. \7\ one (a) corresponding to a 
Naimark-Sacker bifurcation and the other (d) to a pe- 
riod doubling bifurcation. 

As it can be seen from the figures, the vectors (Sg*) 
have several null components. The structural profiles re- 
veal, apparently irrespective of the type of instability, 
that the set S of nodes in the island which are alien 
to instability (white regions), that is, the set of those 
nodes k such that (Sgk) = 0, is a non-zero measure set; 
it is sometimes even larger than the complementary set 
(green area) U = I — S of participating nodes on the un- 
stable eigenvector evolution during a period. We observe 
here that the fragmentation tendency (see discussion on 
islands of disconnected dynamics made above) operates 
also at the level of the tangent space, in the sense that a 
binary partition of the island nodes is well defined at the 
bifurcation (critical) point. Namely, the instability in- 
troduces a partition of the island I = U (B S into (a) the 
set U of nodes that do participate in the instability evo- 
lution in the linear regime, and (b) the complementary 
set iS, of nodes such that single-node perturbations are 
orthogonal to the unstable linear manifold. This drastic, 
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generic fragmentation of the island of periodic activity 
at the linear description level of the bifurcation, is also 
clearly the consequence of the AI competition on the net- 
work of interactions, and we have not seen any deviation 
from this observation in the computations performed (of 
which only two cases are illustrated). In summary, one 
could say that inside the dynamical islands there are com- 
pact substructures (and not single nodes) governing the 
dynamical changes of the whole cluster of nodes. 

The behavior described above suggest the following nu- 
merical experiment: we have explored the responses of 
the different nodes to an external perturbation when the 
system is in a periodic state near a bifurcation point. In 
particular, we force a node by adding an aditional term of 
the form A ■ sin(a;£) (with u> = 2it /T where T the period 
of the unperturbed system) to its equation of motion Q • 
Then we compute, as a function of the forcing amplitude 
A, the evolution of the Floquet eigenvalue /x* responsi- 
ble for the forthcoming bifurcation in the unperturbed 
system. The effects of such a perturbation strongly de- 
pend on whether the perturbed node belongs to the sub- 
set of those identified as leaders [i.e the ones with non 

null component in (6g )). The results obtained for the 
two dynamical islands aforediscussed are shown in figures 
Et c) and I12f f). When the nodes inside the green area 
are perturbed the Floquet eigenvalue fi* significantly de- 
viate (either increase or decrease, we have not been able 
to elucidate when a given change is expected) from the 
values of the unperturbed system. On the other hand, 
the perturbation of the nodes located outside the green 
region does not imply any change to linear stability of the 
whole system. These results illustrate the relevant role 
played by the substructures found by the computation of 

IV. CONCLUDING REMARKS. 

In this paper, we have analyzed the interplay between 
complex topologies and activatory-inhibitory interactions 
driven by a saturated response dynamics of the Michaelis- 
Menten type. The dynamics of the system is very rich 
and exhibits steady, periodic and chaotic regimes that 
in turn lead to the fragmentation of the original sub- 
strate network into a smaller cluster of dynamically ac- 
tive nodes. We have fully characterized these states by 
means of the Lyapunov exponent and the Floquet anal- 
yses as well as the topological features of active islands. 
The reach behavioral repertoire observed is thus a con- 
sequence of the entangled complexity of the system tem- 
poral behavior and the heterogeneous structure of the 
underlying network. 

The emerging dynamics characterized in this work 
could plausible describe at least two relevant scenarios in 
biological systems. On one hand, the dynamics expressed 
in Eq. has been proposed as a way to characterize 
theoretically the individual dynamics of gene expression 
|l3j. In fact, Eq. |Q is the generalization of the suc- 



cessful Random Boolean models [ll|, [12| widely used to 
model gene expression. In this context, two nodes at the 
ends of a link are considered to be transcriptional units 
which include a regulatory gene. One of these end-nodes 
can be thought of as being the source of an interaction 
(the output of a transcriptional unit). The second node 
represents the target binding site and at the same time 
the input of a second transcriptional unit. By studying 
simplified models as the one implemented here — the in- 
trinsic complexity of the problem does not allow for a 
complete and detailed description of real gene dynamics 
— , one can infer the region of the parameter space (i.e. 
(p, h)) that better describes gene networks. The latter 
seems possible due to latest developments in microarray 
technologies, biocomputational tools, and data collection 
software. 

A second scenario where the results obtained apply is 
reaction kinetics in metabolic networks. In metabolic 
systems, a very rich behavioral repertoire is well doc- 
umented [23, as for instance, the oscillations observed 
in the concentration of certain chemicals in biochemical 
reactions such as glycolysis. The system of differential 
equations, Eqs. represents one of the most basic bio- 
chemical reactions, where substrates and enzymes are in- 
volved in a reaction that produces a given product. In 
this context, there are several important issues as how 
fast the equilibrium is reached, how the concentration 
of substrates and enzymes compare, etc. Besides, it is 
known that in a large number of situations, some of the 
enzymes involved show periodic increments in their ac- 
tivity during division, and these reflect periodic changes 
in the rate of enzyme synthesis. This is achieved by reg- 
ulatory mechanisms that necessarily require some kind 
of feedback control as that emerging in our model. The 
interesting point here is that the real topological features 
of the underlying metabolic network 23] have not been 
taken into account in studies performed so far. As this 
work shows, they have important bearings in the corre- 
lation between structure and the observed dynamics. 

Finally, on more general theoretical grounds, we antic- 
ipate several features of interest such as the fragmenta- 
tion of the original network according to the dynamical 
states of the nodes, multistability and different routes to 
chaotic behavior within the same system. The first of 
these points is particularly relevant since it may indicate 
that in networks of dynamical units, the topology ob- 
served can be the result of a given network state hiding 
a larger substrate whose topological properties are com- 
pletely different at a local level. Of particular interest 
is also the result gathered in the last part of the work, 
namely, the existence of an additional substructure inside 
dynamical islands determined by the different responses 
of nodes to external perturbations. This points to the 
central issue in many biological processes of what subset 
of nodes are the most important in order to sustain (or 
break) the system's robust functioning. In summary, the 
characterization of models where nonlinearity and spatial 
complexity coexist yields new results missed when only 
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one of these ingredients is present and opens the path to 
a better comprehension of biological processes and the 
dynamics of networks of nonlinear dynamical units. 
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